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Abstract 

In this paper we study a Hamiltonization procedure for me- 
chanical systems with velocity-depending (nonholonomic) 
constraints. We first rewrite the nonholonomic equations 
of motion as Euler-Lagrange equations, with a Lagrangian 
that follows from rephrasing the issue in terms of the inverse 
problem of Lagrangian mechanics. Second, the Legendre 
transformation transforms the Lagrangian in the sought-for 
Hamiltonian. As an application, we compare some varia- 
tional integrators for the new Lagrangians with some known 
nonholonomic integrators. 

Keywords: nonholonomic systems, Lagrangian and Hamil- 
tonian formalism, inverse problem, geometric integration. 

1 Introduction 

Many interesting mechanical systems are subject to addi- 
tional velocity-dependent (i.e. nonholonomic) constraints. 
Typical engineering problems that involve such con- 
straints arise for example in robotics, where the wheels 
of a mobile robot are often required to roll without slip- 
ping, or where one is interested in guiding the motion of 
a cutting tool. 

The direct motivation for our paper ||2l was to be found 
in interesting results that appeared in Q, where the au- 
thors propose a way to quantize some of the well-known 



classical examples of nonholonomic systems. On the way 
to quantization, the authors propose an alternative Hamil- 
tonian representation for those nonholonomic systems. 
However, the "Hamiltonization" method introduced in O 
can only be applied to systems for which the solutions are 
already known explicitly. 

Nonholonomic systems have a more natural descrip- 
tion in the Lagrangian framework. In |2|, we explained 
how one can associate to the nonholonomic equations of 
motion a family of systems of second-order ordinary dif- 
ferential equations and we applied the conditions of the 
inverse problem of the calculus of variations on those as- 
sociated systems to search for the existence of a regular 
Lagrangian. (The inverse problem of the calculus of vari- 
ations deals with the question of whether or not a given 
system of second-order differential equations is equiva- 
lent with the Euler-Lagrange equations of a yet to be de- 
termined regular Lagrangian, see e.g. |T2l). If such an un- 
constrained regular Lagrangian exists for one of the asso- 
ciated systems, we can always find an associated Hamilto- 
nian by means of the Legendre transformation. Since our 
method only made use of the equations of motion of the 
system it did not depend on the knowledge of its explicit 
solutions. 

A system for which no exact solutions are known can 
only be integrated by means of numerical methods. In 
addition to the above mentioned appUcation to quanti- 



zation, our Hamiltonization method may also be useful 
from this point of view. Numerical integrators that pre- 
serve the underlying geometric structure of a system are 
called geometric integrators. A geometric integrator of a 
Lagrangian system uses a discrete Lagrangian that resem- 
bles as much as possible the continuous Lagrangian (see 
e.g. mi). On the other hand, the succes of a so-caUed 
nonholomic integrator (see e.g. p] W\) relies not only on 
the choice of a discrete Lagrangian but also on the choice 
of a discrete version of the constraint manifold. It seems 
therefore reasonable that if a free Lagrangian for the non- 
holonomic system exists, the Lagrangian integrator may 
perform better than a nonholonomic integrator with badly 
chosen discrete constraints. 

In the next section we recall the set-up and the main 
results of our paper |2|. In section [3] we compare some 
nonholonomic and variational geometric integrators for a 
few of the classical nonholonomic systems. In section]?] 
we indicate some ideas on how we wish to extend the re- 
sults of this paper 



2 A class of nonholonomic systems 

We will consider only a certain class of nonholonomic 
systems on M": We will assume that the configuration 
space of the system is a space with coordinates {ri,r2,Sa), 
that the Lagrangian of the system is given by the function 



(1) 



and that the nonholonomic constraints are all of the form 



(2) 



The nonholonomic equations of motion follow from 
d'Alembert's priciple (see e.g. [IJ). For systems in our 
class they are given by the equations 
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together with the constraint equations (|2]l. After eliminat- 
ing the Lagrange multipliers by means of the constraints, 
one gets 



r2 - 



(3) 



where N{r\ ) — (/2 + La^a^a) - is related to the invariant 
measure of the system and K = Lp^p^P'^p. with Ap = 

Some basic examples of nonholonomic systems that lie 
in this class are the following ones. The classic exam- 
ple of a nonholonomically constrained free particle has a 
Lagrangian and constraint given by 



1 



and z+xy~Q. 



A knife edge on a horizontal plane corresponds physically 
to a blade with mass m moving in the xy plane at an angle 
(|) to the X-axis. Its Lagrangian and constraint are given by 



1 



1 



L^-m{x^+f) + -J^^ 



and 



- jcos( 



0. 



Also the vertically rolling disk is an example in our class. 
The assumption that the disk rolls without slipping over 
the plane gives rise to nonholonomic constraints. Let R be 
the radius of the disk. If the triple [x^y^z — R) stands for 
the coordinates of its centre of mass, (p for its angle with 
the (jCjz) -plane and for the angle of a fixed line on the 
disk and a vertical line, then the nonholomic constraints 
are of the form 

x = /?cos(p0 and y = Rsin(pQ. 

The Lagrangian of the disk is 



M{x^+f 



1 



where / — jMR^ and J = ^MR^ are the moments of iner- 
tia and M is the total mass of the disk. For the vertically 
rolling disk is a constant and K — 0. 



z 




Figure 1 : The vertically rolling disk 



Finally, also the examples of the mobile robot with fixed 
orientation and the two- wheeled carriage (see e.g. |8 1) lie 
within our class. 

The equations of motion Q are a mixed set of first- and 
second-order differential equations. On the other hand, 
the Euler-Lagrange equations 

dt \dc{) ^ 

of a regular Lagrangian L are second-order differential 
equations (only) [The tilde in L will always denote that the 
Lagrangian is free, and that it should not be confused with 
the original Lagrangian L of the nonholonomic system.]. 
We therefore need a way to associate a second-order sys- 
tem to our nonholonomic system. One possible choice of 
doing so is the system 

r n=0, 

J h = -N^Krir2, (4) 
[ s^ = -{A'^-N'-KAa)hr2- 

The above second-order system has the property that its 
solution set contains, among other, also the solutions of 
the nonholonomic dynamics ([3]l when restricted to the 
constraints. Another choice for an 'associated system' 



with the same property is e.g. 

'n=0, 

Sa = {A'^-N^KAa)h 

(no sum over a). It is clear, that there are in fact an infinite 
number of such associated second-order systems, but we 
will concentrate in this paper on the above two. For some 
other possible choices, see 

Proposition 1. L There does not exist a regular La- 
grangian whose Euler-Lagrange equations are equivalent 
with the second-order system Q (for the classical exam- 
ples cited above). 

2. The Euler-Lagrange equations of the Lagrangian 

~^-W^ + ^{c2-+l.C^^\ (C„^0) (6) 

are equivalent with the second-order system (|5]l. If the 
invariant measure density is a constant, then also 

L=^/i.?4/..f+^i:.p^ (c„^o) (7) 

is a regular Lagrangian for the system (j5]l. 

Proof. We give here only an outline of the method we've 
used to prove the statements. For full details, see fT\. As- 
sume we are given a system of second-order ordinary dif- 
ferential equations 

The search for a regular Lagrangian is known in the litera- 
ture as 'the inverse problem of the calculus of variations', 
and has a long history (for a recent survey on this history, 
see e.g. |10| and the long list of references therein). In 
order for a regular Lagrangian L{q, q) to exist we must be 
able to find functions gij{q,q), so-called multipliers, such 
that 



It can be shown |l6][T2l that the muhipHers must satisfy 

^8U _ ^gik . 



where 



and 



detfe;)7^0, gji^gij, 



dq^ dqj 



The symbol F stands for the vector field q'd^i + f'd^i that 
can naturally be associated to the system q' = f'{q,q). 
Conversely, if one can find functions gij satisfying these 
conditions then the equations q' ~ f are derivable from a 
regular Lagrangian. Moreover, if a regular Lagrangian L 

can be found, then its Hessian t^^^t — r is a multiplier. 

dq'dqJ 

The above conditions are generally referred to as the 
Heknholtz conditions. They are a mixed set of coupled 
algebraic and PDE conditions in {gij). We will refer to 
the penultimate condition as the 'V- condition,' and to the 
last one as the 'O-condition.' The algebraic <I>-conditions 
are of course the most interesting to start from. In fact, we 
can easily derive more algebraic conditions (see e.g. |5 1). 
For example, by taking a F-derivative of the O-condition, 
and by replacing F(^,y) everywhere by means of the V- 
condition, we arrive at a new algebraic condition of the 
form 

where 

(V4>)} = F(<i>^)-v;,<i>7-v74>;„. 

This (V<I>)-condition will, of course, only give new infor- 
mation as long as it is independent from the <I>-condition 
(this will not be the case, for example, if the commutator 
of matrices [<I>, V<I>] vanishes). One can repeat the above 
process on the ( V«l>)-condition, and so on to obtain possi- 
bly independent (V . . . VO) -conditions. A second route to 
additional algebraic conditions arises from the derivatives 
of the <I>-equation in ^'-directions. One can sum up those 
derived relations in such a way that the terms in 3(^^17 



disappear on account of the symmetry in all their indices. 
The new algebraic relation in gij is then of the form 

where 

As before, this process can be continued to obtain more 
algebraic conditions. Also, any mixture of the above men- 
tioned two processes leads to possibly new and indepen- 
dent algebraic conditions. Once we have used up all the 
information that we can obtain from this infinite series of 
algebraic conditions, we can start looking at the partial 
differential equations in the V-conditions. 

Let us now come back to the second-order systems (j4]) 
and (j5]l at hand. The proof of the proposition relies on the 
fact that for the first systems (j4]i, the only matrices (gij) 
that satisfy the first few algebraic conditions must be non- 
singular On the other hand, the two Lagrangians for the 
system (j5]l follow from an analysis of the Helmholtz con- 
ditions with carefully chosen anszatzes. For more details, 
see 0. □ 

Remark that the Lagrangians are not defined for ri = 0, 
and we will in general exclude the solutions with that 
property from the further considerations in this paper. 
Any regular Lagrangian system with Lagrangian L can be 
transformed into a Hamiltonian one, by making use of the 
Legendre transformation 

{q\q') {q\pi = ^). 

The corresponding Hamiltonian is then 

H = Piq'-L. 

Similarly, the Legendre transformation maps the con- 
straints, viewed as a submanifold of the tangent manifold, 
onto a submanifold in the cotangent manifold. 

Proposition 2. Using the Legendre transformation, the 
Hamiltonian that corresponds to the Lagrangian (j6]) is 
given by 



The corresponding constraints are 

CiPa = -CaPl- 

If is constant, the Hamiltonian that corresponds to the 
Lagrangian (jTjl is 



H 




and the constraints transform into 

hNriPa + aaP2=Q, 
where ri = {pi + \NY.aAaPl./aa) Ih- 

In |2 9 1 we explain how the above Hamiltonians can be 
directly derived from Pontryagin's Maximum principle. 

3 Geometric integrators 
3.1 Set-up 

As we explained in the introduction, there are now two 
ways to compute a numeric approximation of a solution of 
a system in our class: we can use either a nonholonomic 
integrator for the original Lagrangian ([T]i and constraints 
(j2]i , or we can use a variational integrator for one of the 
Lagrangians (|6]l and (j7]| we have found in Proposition [T] 
Let us come to some details. 

Geometric integrators are integrators that preserve the 
underlying structure of the system. In particular, varia- 
tional integrators are integrators that are derived from a 
discrete version of Hamilton's principle. From this dis- 
crete variational principle one obtains the so-called dis- 
crete Euler-Lagrange equations as follows. For a me- 
chanical system with Lagrangian L, one needs to choose 
a discrete Lagrangian Ld{qi,q2) (a function on Qx Q 
which resembles as close as possible the continuous La- 
grangian). A solution q{t) is then discretised by an array 
qii which are the solutions of the so-called discrete Euler- 
Lagrange equations 



DiLd{qk,qk+i)+D2Ld{qk-i,qk) = 0. 



(8) 



These integrators preserve the symplectic and conserva- 
tive nature of the algorithms. It is important to realize that 



a different choice for the discrete Lagrangian may lead 
to a different geometric integrator The presence of ad- 
ditional holonomic constraints (i.e. 'integrable' nonholo- 
nomic constraints) can be included by introducing La- 
grange's multipliers. 

On the other hand, for a nonholonomic integrator of 
a nonholonomic system with Lagrangian L and con- 
straints (0"{q)q" — 0, we need to choose both a discrete 
Lagrangian Lj and discrete constraint functions co^J on 
Qx Q. The nonholonomic discrete equations are then 

J D\Lci[qk,qk+\) +D2Ld{qk-i,qk) = {h)a(£i"{qk), 
\ (O"/qk,qk+l)=0. 

(9) 

Usually, if g is a vector space, one takes the discretization 
in one of the following ways (for certain a and certain h): 



Ld{qi,q2) 



«d(?l>?2) = 



L[q={l - a)q\ + aq2,q - 



lb) 



12 -qi 



<(^ = (l-a)^i+a^2)^S:^. (11) 



For the rest of the paper, we will concentrate on this dis- 
cretization procedure. There are, however, many more 
possibilities to obtain a discrete Lagrangian and discrete 
constraints. For example, one could take a symmetrized 
version of the above procedure and use discrete La- 
grangians and discrete constraints of the form 

Ldiqi,q2) = (^q = {l-a)qi+aq2,q = ^ ^' 

1 / , , q2^q\ 

+ -L\q = aqi + {l-a)q2,q= 



oi"i{quq2) = ^< (g = (1 - + «g2) ^ ^' 

+ (q = aqi + (1 - a)q2) ■ 

Also, if the system is invariant under a symmetry group, it 
is advantageous to construct the integrator in such a way 
that the discrete system inherits as many as possible of 
those symmetry properties, see e.g. [4 |. 

The bottom line of the next sections is the following 
one. If a free Lagrangian for the nonholonomic system 
exists, it seems reasonable that the Lagrangian integrator 



may perform better than a nonholonomic integrator with 
badly chosen discrete constraints. In the next sections, we 
will test this conjecture on a few of the classical exam- 
ples in our class: the vertically rolling disk, the knife edge 
and the nonholonomic particle. It will be convenient that 
for those systems an exact solution of the nonholonomic 
equations ([3]) is readily available. 

3.2 The vertically rolling disk 

For the vertically rolling disk, we have (ri,r2,Sa) = 
(9, Q,x,y). It is well-known that the solutions of the non- 
holonomic equations with initial conditions = (p(0) 7^ 
and Me = 0(0) are all circles with radius R{uq/u^): 



e(f) = Mef + Qo 
Me 



x(f) 



yit) = 



Me 



(p(f) = M,pf + (po, 

/;sin((p(f)) +X0, 
/;cos((p(f))+yo- 



(12) 



Let us put for convenience M = I and R = I and there- 
fore / = 2 and J = \- With that the (nonholonomic) La- 
grangian and constraints are simply 



1 



-cp , i = cos90, 3> = sin(p9. 



We will first compute the solution of the discrete nonholo- 
nomic equations (j9|l with the discrete Lagrangian ( 10 1 and 
the discrete constraints 



11 



Second, since the vertically 
rolling disk is one of those examples with a constant in- 
variant measure density A^, we can choose a Lagrangian 
from the second type (|7]). The simplest choice is probably 



L= 1/2 [ (p2+e2- 



•7 -2 \ 

+ y ] 

cos((p)(p sin((p)(p / 



(13) 



We now investigate the variational integrator of this La- 
grangian, where the discrete Lagrangian is given by ( [T0| . 
We will fix h (changing it did not have a significant effect) 
and only concentrate on what happens if we keep a vari- 
able. In figure|2]we have plotted the situation for a = 0. 
For a given set of initial positions (jccycQojtpOjQij^Pi) 
the other initial conditions were chosen in such a way that 




Figure 2: Vertically rolling disk with a = 0. 



the solution lies initially on the discrete constraint mani- 
fold, i.e. in such a way that 

xi ::=xo+cos(po(6i -9o), yi =yo + sin(po(ei -60). 

Unlike the nonholonomic integrator (in grey with circle 
symbols) the variational integrator (in black with cross 
symbols) does not show a strong spiral-type solution but a 
circular path. It is true, however, that the variational solu- 
tion deviates from the circle predicted by the initial con- 
ditions of the solution (12 1 (in grey in figure |2|i. However, 
since any circle is determined by 3 of its points, we can 
find a better match for the circle the variational discrete 
solution follows by considering the outcome {xkj,yki) at 
three different times and by solving the three equations 



{xk_-Af + ixk^~Bf^C^ 

for {A,B,C). If we do so, we obtain the matching circle 
(in dots) in figure |3] 

It is well-known that the energy 

E^\{x'+f) + \B' + l<^\ 

is conserved along the solutions ( [T2] i of the nonholonomic 
equations of motion. In figure|4]we investigate the perfor- 
mance of the two integrators on the energy function. The 
discrete version of the energy is the function we get by 
substituting, as usual, 1 / h{qi; — qk- 1 ) for q in the function 
E above. The straight line in figure]?] is the energy level 





Figure 3: Vertically rolling disk: circular path. 



fied' variational integrator which does conserve the con- 
straints. This integrator considers the constraints as a con- 
stant along the (nonholonomic) motion. That is, it will 
use the variational discrete Lagrange equations (jSjl for the 
variables 9 and (p (for the free Lagrangian L given in ( 13 1), 
but not the corresponding equations for x and y. To get a 
full system of equations, we supplemented this with the 
discrete constraints (i)"^{qkTqk+\) — which can be writ- 
ten in terms of Xjt+i and jjt+i . Figure[5]shows the modified 
integrator for a = (with box symbols). The circle in that 
figure is the one we had before, i.e. the one that matches 
the variational integrator. It shows that the modified inte- 
grator has the same circular behaviour as the variational 
integrator, and on top, it keeps the constraints conserved, 
see the box symbols on the zero level in figure [6] 



predicted by the initial conditions. It is clear that the vari- 
ational integrator (with crosses) does a better job than the 
nonholonomic one (with circles). 




Figure 4: Vertically rolling disk: the energy. 

By construction the nonholonomic integrator conserves 
the constraints and the variational integrator does not. 
Indeed, in figure |6] we have plotted the constraint i — 
7?cos((p)9 = 0. Positive is that, although the variational 
integrator does not conserve this constraint, it reasonably 
oscillates around the zero level. Moreover, there is a 
method to fix this problem. We can introduce a 'modi- 



k 

i 



Figure 5: Vertically rolling disk: the modified integrator. 

Finally, figure |7] shows the effect of changing the pa- 
rameter a. The results for the variational integrator (in 
black with cross symbols) remain accurate and more or 
less unchanged. For the nonholonomic integrator (in grey 
with circle symbols) the effect of changing a is that the 
inward spiral becomes an outward spiral. At some point 
(here a = 1/2) the variational and nonholonomic integra- 
tor have the same accuracy. 

3.3 The knife edge 

As was the case with the vertically rolling disk, also the 
solutions of the knife edge form a circular path in the 
(x,3')-plane. Continuing the analogue with the previous 




Figure 6: Vertically rolling disk: the constraints. 



Figure 8: The knife edge with a = 





Figure 7: Vertically rolling disk with a = 1/3,1/2,1, 
repectively. 

example, the nonholonomic integration (in grey with cir- 
cle symbols) results in a spiral, while the variational in- 
tegration (in black with crosses) follows more closely the 
circular path, see figure |8] 

3.4 The nonholonomic particle 

The function 



- 1 , 
2 



z 

XX 



is a free Lagrangian for the nonholonomic particle. In 
each of the figures|9]and 10 the dashed black curve repre- 
sents the exact solution, the thick black the variational so- 
lution and the thick grey the nonholonomic solution. The 
figures show that both the variational method and the non- 
holonomic one do not give very accurate solutions. How- 



ever, changing the parameter a does not seem to affect the 
variational solution as much as it does the nonholonomic 
one. Indeed, the variational solution remains more or less 
of the same accuracy for the different a-values. On the 
other hand, the nonholonomic solution can be made more 
or less accurate by changing a. It seems that the best 
accuracy is reached somewhere in the neighbourhood of 
a = 1/3, but how could one have guessed this before- 
hand? Remark also that this value is not same as the the 
best choice we had found for the nonholonomic integrator 
of the vertically rolling disk (where a = 1/2 gave the best 
accuracy). 

3.5 Preliminary conclusion 

In each of the discussed examples the variational integra- 
tor (with one of the Lagrangians ^ and (|7]i of proposi- 
tion [TJ seemed to give better results than the known non- 
holonomic integrators. Unlike the outcome for the non- 
holonomic integrator, the results for the variational inte- 
grator seemed to be independent of or, at least, stable un- 
der changing the parameter a. Needless to say, the results 
above are, of course, very partial and are they are only in- 
tended to motivate further investigation on this topic. For 
example, we need to check if more involved discretiza- 
tion procedures, such as the ones mentioned at the end of 
section[3j\ demonstrate the same behaviour as the one we 
have encountered so far. 





mobile robot with a fixed orientation 

L = \m{x^ +f) + i/e^ + ^7xj;2 „ iOsin\|/, 
x — RcosQ^, y = RsinQ^, 

(an example that also appears in the paper ||4l) or the knife 
edge on an inclined plane, where 

L= -m{x^ +y^) + -J^^ +mgxsina, isin(|) = jcoscj). 

In more general terms, such systems have a Lagrangian of 
the form 

1 




L^-{hri+l2ri+bs')-Vir2) 
and a constraint of the form 

s = -A(ri)r2, 

\\and we can we can consider associated second-order 
slpquations, in a way that is analogous to the way we ar- 
ived at the second system Q before: They are now of 
e form 



Figure 9: The nonholonomic particle: xy-, xz- and yz- 
solution with a = 0. 



4 Further systems 



The class of nonholonomic systems treated above is very 
restricted. The reason is, of course, that the search for a 
solution of the inverse problem of the calculus of varia- 
tions (in the proof of proposition [T]i is too hard and too 
technical to be treated in the full generality of a nonholo- 
nomic systems with an arbitrary given Lagrangian and ar- 
bitrary given constraints. Also, since there are infinitely 
many possible choices for the associated systems, it is not 
clear from the outset which one of them will be varia- 
tional, if any. 



n =0, 

h = r2(ri)nr2 + f2(n,'"2), 



(14) 



Remark that compared to the equations (j5]l, the presence 
of the extra potential brings the terms f,(ri,r2) into the 
picture. A first result is the following. 

Proposition 3. There does not exists a regular Lagrangian 



for the second order systems (14 1 



Proof. As before, the proof follows from a careful analy- 
sis of the algebraic conditions which can be derived from 
the Helmholtz conditions. □ 

For systems with more than one constraint, the result is 
still open. Remark that the proposition does not exclude 
the existence of an other variational 'associated' system. 



For these reasons, future extensions of the obtained re- 
sults will strongly depend on well-chosen particular new AcknOWledgllieiltS 

examples. For example, we could try to find a free La- 
grangian for a nonholonomic system with a potential of 
the form V{r2). Typical examples of such systems are the 
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Figure 10: The nonholonomic particle: xy-solution with 
a = 0, 1/5, 1/3, 1/2,2/3,4/5, respectively. 



